Finite Size Fluctuations in Interacting Particle Systems 
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Fluctuations may govern the fate of an interacting particle system even on the mean-field level. 
This is demonstrated via a three species cyclic trapping reaction with a large, yet finite number 
of particles, where the final number of particles Nj scales logarithmically with the system size A*', 
Nf ~ InA'^. Statistical fluctuations, that become significant as the number of particles diminishes, 
are responsible for this behavior. This phenomenon underlies a broad range of interacting particle 
systems including in particular multispecies annihilation processes. 
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I. INTRODUCTION 

Balls-in-boxes (urn) models provide a handy labora- 
tory for studying conceptual issues. The celebrated Eren- 
fest model for instance, has led to the reconcil- 

iation of the reversibility and the recurrence paradoxes 
with Boltzmann's H theorem. Recent examples include 
urn models of aging (J, |M 1^ and of discretized quantum 
gravity . Urn models have been studied extensively in 
probability theory loin^ and have found applications 
ranging from biology |l ij to computer science 0, 0| . 

Perhaps the most well-known manifestation of the role 
of fluctuations in stochastic processes is the Eggenberger- 
Polya urn model One starts with two marbles of 

different colors, draws a marble randomly and puts it 
back together with another marble of the same color. 
When the process is repeated indefinitely, the fraction 
of marbles of a given color saturates at some limiting 
value. The corresponding limit is a random variable that 
is uniformly distributed between zero and one. Thus, 
initial fluctuations are locked in, a striking example of 
the lack of self-averaging. 

Inspired by this example, our goal is to quantify fl- 
nitc size fluctuations in interacting particle systems us- 
ing urn models. Finite-size corrections are important 
because they govern for example how a system con- 
verges to the thermodynamic limit, yet they remain 
largely unresolved even in elementary stochastic pro- 
cesses 111 111 in 111 IIS Ei. 



We study the role of fluctuations using a three color 
cyclic urn model. Initially, the urn contains three dif- 
ferent types of balls. Then, two balls are picked ran- 
domly. If they are dissimilar, following a cyclic rule, one 
of the balls is returned to the urn and the other is re- 
moved from the system. This urn model is different from 
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the Eggenberger-Polya type models in two ways. First, 
the number of balls is decreasing rather than increasing. 
Second, our model is nonlinear because picking two balls 
rather than one implies that different type balls interact 
with each other. 

We deflne the state of the system when one of the 
species is depleted to be the flnal state. Starting from 
the natural initial conditions where there are N balls of 
each type, we ask: "How many balls are there in the flnal 
state?". Our main result is that the average number of 
balls in the final state Nf scales logarithmically with the 
system size N . 

Statistical fluctuations are ultimately responsible for 
this behavior. As long as the system contains enough par- 
ticles (balls), the average number of particles faithfully 
characterizes the state of system. However, as particles 
deplete, the uncertainty with respect to how many par- 
ticles remain grows and moreover, it governs how many 
particles are flnally left. 

This phenomenon and the mechanism underlying it are 
generic to interacting particle systems with a decreasing 
number of particles. We demonstrate this by consider- 
ing multi-species annihilation processes with q species. 
In the three species model, there is again a logarithmic 
enhancement of the variance in the number of particles 
over the average number of particles, thereby leading to 
a logarithmic growth law. In general, there is an al- 
gebraic enhancement as long as the number of species 
is small enough, and consequently, an algebraic growth 
law. Otherwise, when the number of different species is 
large enough, statistical fluctuations are negligible and 
the number of remaining particles is of order one. In the 
most general case when p balls are drawn from the urn, 
the critical number of species is qc = 2p — 1. 

The rest of the paper is organized as follows. In Sec. II, 
we give a nontechnical presentation of the cyclic trapping 
model and highlight the basic result. We then analyze 
the model in detail and obtain the number density fluc- 
tuations by employing the 1/N expansion (Sec. HI). The 
q-species annihilation model is treated in Sec. IV. We 
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conclude with a few open questions (Sec. V). 



II. CYCLIC TRAPPING REACTIONS 

Let us first define the model. Initially, the urn contains 
balls of three different types, denoted by A, B, and C. The 
balls interact via a cyclic trapping reaction. Two balls 
are taken randomly out of the box. If they are different, 
then one of the balls is returned to the urn according to 
the cyclic rule 

A + B^B, B + C^C, C + A^A, (I) 

and the second ball is discarded. This elemental step 
is repeated until one of the species becomes extinct. 
This reaction scheme is reminiscent of the Lotka-Volterra 
cyclic food chain (or rock-paper-scissors) model, widely 
used in ecology and game theory [Tlll2lll2ll2ll2i| . 

The state of the system is fully specified by the num- 
ber of particles of each type in the urn: n, m, and /, 
corresponding to particles of type A, B, and C. The dy- 
namics is clearly mean-field: every dissimilar pair of par- 
ticles is equally likely to interact. Therefore, the tran- 
sition {n^m,l) — !■ (n — l,m, Z) occurs with probability 
nm/ (nm + ml + In) and similarly for the two other tran- 
sitions. 

We start with the initial condition n = m = I ^ N 
and stop the process when one species becomes extinct. 
Surprisingly, the number of balls in the final state scales 
logarithmically with the system size: 



N 



f 



IniV. 



(2) 



Numerical simulations are consistent with this behavior 
(Fig. 1). We verified numerically that the scale In fully 
characterizes the final number of particles. The distribu- 
tion of the final number of particles approaches a (non- 
trivial) limiting distribution when the final number of 
particles is normalized by huN. Thus, the final number 
of particles in a non-self-averaging quantity. 

The problem is essentially combinatorial, and in prin- 
ciple, it can be addressed by weighing all possible his- 
tories with the appropriate transition probabilities |25j |. 
It proves fruitful, however, to treat the process dy- 
namically. Choosing the rate nra/N for the transition 
(n, m, /) ^ (n, m — 1, 1) is consistent with the above mi- 
croscopic rules. Moreover, it leads to A-independent dy- 
namics in the thermodynamic limit. 

The number densities a — {n)/N, b = {m)/M, and 
c = {l)/N evolve according to the rate equations 



da 
'dt 



-ab. 



db 
dt 



= -be 



dc 
di 



(3) 



Since initially a(0) = 6(0) = c(0) = 1, the number den- 
sities remain equal throughout the entire process a{t) ~ 
b{t) = c{t) = p{t) with 




FIG. 1: The average total number of particles in the final 
state as a function of the system size. The data represents an 
average over 10^ realizations of the cyclic trapping reaction 
process Q. 



Naively assuming that throughout the process, fluctua- 
tions in the number density are much smaller than the 
mean leads to the conclusion that the final number of 
particles is of the order unity, Nf ~ 0(1). The corre- 
sponding terminal time scales linearly with the system 
size, tf N . Below, we show that this assumption does 
not hold when the total number of particles becomes suf- 
ficiently small. 

The logarithmic growth in the number of particles can 
be deduced from the fluctuations in the number density. 
In the thermodynamic limit, we expect that to leading 
order in A, both the total number of particles and the 
variance in the number of particles are proportional to 
the system size 



Ap, 
Acr^ 



(5) 



We term the intrinsic variance. In Sec. IIL we shall 
utilize the van Kampen 1/A expansion |26l l27l| to show 
that asymptotically, the ratio between the intrinsic vari- 
ance and the density grows logarithmically with time 



Ini. 
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(6) 



Thus, fluctuations eventually become larger then the 
density. Of course, when they are comparable with 
the density, extinction is possible. Hence, the number 
density Q characterizes the particle number only up 
to a time scale tf obtained from the validity criterion 
Np{tf) ^ ^JNaHtf). The terminal time is therefore 



tf - A(lnA)" 



(7) 



p{t) = {i + ty 



(4) 



Using Nf ~ Np{tf) we arrive at our main result 
Note that In A is the leading contribution. The sub- 
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leading contribution ln(lnA'^) corresponding to nested 
logarithms is tacitly ignored. 

III. PARTICLE NUMBER FLUCTUATIONS 

Fluctuations in the particle number are studied by ex- 
panding the master equation in inverse powers of N and 
keeping only the leading order terms (large- expansion) 
[2^ . The probability P{n, m, I; t) that the particle num- 
bers are n, m, and I at time t obeys the master equation 



-P(n, TO, I) = {Cab + Cbc + Cca) P{n, m, (8) 
at 

with the initial condition Po{n,m,l) = 6n,NSm.N^i,N- 
The operator Cab is 

CABP{n,m,l) N-'^{Aa - 1) [nmP{n,m,l)] ; (9) 

the operators Cbc a-nd Cca are defined via similar for- 
mulas. The difference operator A raises the respective 
variable by one, e.g. 

AAf{n,m,l)^ f{n + l,m,l). (10) 

Since in the thermodynamic limit, averages as well as 
variances grow linearly with the system size as in Eq. |(SJ), 



where (c.t.) denotes the two terms obtained by cyclic 
transposition of the displayed term on the right-hand 
side. This master equation contains terms of various or- 
ders in N. The order N^^'^ terms vanish because the 
densities satisfy the rate equations (01 . The next leading 
order term gives the evolution operator 

MABF = pda[{a + l3)F] + ^p^F^^. (13) 

Explicitly, the Fokker-Planck equation (|ll|l is 

Ft = p[da,ia + f3) + dpil3 + j) + d^{j + a)]F 

+ ^p\Fac + F0f3 + F^^). (14) 

This Fokker-Planck equation is subject to the initial con- 
dition Fo{a, (3, 7) = (5(a) 6{(3) S{j). 

Moments of the probability distribution F{a, (3, 7) di- 
rectly follow from H14() . One simply multiplies this 
Fokker-Planck equation by the desired powers of a, (3, 
and 7, and integrates (by parts) with respect to these 
three variables. Due to symmetry, there is essentially 



we introduce the transformation P{n, m, I) — s- F{a, P,"f) 
with 

n^Na + N^^^a, Nb + N^^^P, I = Nc + N^^^j. 

The intensive (random) variables a, /?, and 7 are A''- 
independent. These variables simply characterize fluctu- 
ations in the respective particle numbers. 

To find out how the distribution F(a, /?, 7) evolves with 
time, we write 

Ft = {MAB+MBC+McA)F. (11) 

It suffices to compute the evolution operator A4abF; 
the two other operators are obtained by cyclic trans- 
position. To obtain the evolution operators, we re- 
place the distribution P by F in (|SJ) and convert dif- 
ference equations into differential ones by expanding dif- 
ference operators and keeping up to second order terms, 
e.g., ^ 1 + dn + \dnn- Similarly, we replace deriva- 
tives with respect to n with derivatives with respect to 
a using 9„ = N^^/^da- The time derivative becomes 
dt - N^/^ada - N^/^bdp - N^/^cd^ where the overdot 
denotes differentiation with respect to time. These trans- 
formations lead to 



(12) 

I 

one first moment: (a); two second moments; (a^), (a/3); 
three third moments: (a^), {a^f]), {a/Sj); etc. The first 
moment satisfies -^{o) — —'2p{a) and since it vanishes 
initially, (a) = 0. The two second moments are coupled 

^ = -2p{a^)-2p{aP)+p^, (15) 

^ = -p{a')-3p{aP). 

These equations are inhomogeneous, so despite of the 
vanishing initial conditions (a^) = (a/3) = 0, the solu- 
tions are non-trivial. 

Writing U = {a^} + 2 (a/3) and V ^ (a^) - (a/3), we 
separate the above equations 

— = -4pU + p^ (16) 

-M - + ^ ■ 
Using the number density Q), we obtain the explicit ex- 



I 

Ft - N^/^{dFa + bFp + cF^) = {N'^da + ^N-^/^dca)[iNa + N^/^a){Nb + N'^''^(3)F] + (c.t.) 
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pressions 



U 



[{i + tr'-ii + ty 



(17) 



V = (l + t)-Mn(l+t). 



Physically, U = {a{a + /3 + 7)) quantifies the correla- 
tion between the single particle number n and the total 
particle number n + m + I, while V = {a{a — /3)) quan- 
tifies the correlation between the particle number n and 
the number difference n — m. Intuitively, we expect that 
the quantity V is larger than U. For a sufficiently large 
system, it may be arbitrarily larger. 

One of the two second moments is the intrinsic variance 



explicitly. 
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P 



ln(l 



X 1 1, 



(18) 



The other (normalized by the density) second moment 
quantifies cross-correlations between different species 
numbers 



(aP) 



ln(l 



' 3 3^ 



(19) 



The quantity (a/3) is always negative and therefore, fluc- 
tuations between different particle numbers are anti- 
correlated. Asymptotically, {a(3) ~ —0^/2 with 



2 

3^" 



Int. 



(20) 



Another important consequence of the structure of the 
Fokker-Planck equation is that the multivariate distribu- 
tion P(n, TO, I) is Gaussian and fully characterized by the 
first and second order moments. This is the case because 
the first order derivatives in (|14() have linear coefficients 
|26l |. As a result, the individual particle number distri- 
bution is also Gaussian 



1 



exp 



{n - Npf 



2Na^ 



(21) 



IV. MULTISPECIES ANNIHILATION 

We have examined the question "how many particles 
remain in the final state?" in a number of other interact- 
ing particle systems where depletion or extinction occurs 
and found that generically, fluctuations play an impor- 
tant role. Using the same validity criterion, and utilizing 
the van Kampen's 1/A'' expansion, one can determine the 
final number of remaining particles. 

We demonstrate this for multi-species annihilation pro- 
cesses. Initially, the urn contains q types of balls and the 
initial number of each species is equal to N . For instance, 
when q — 



A + S -> 0, 



C + A^O. 



(22) 



This process, introduced by ben-Avraham and Redner, 
was studied primarily in low spatial dimensions via a 
number of numerical and analytical techniques, yet it is 
still not fully understood HHHl^. 

The parameter q is in principle integer. However, it is 
still sensible to treat it as a continuous variable in the 
range 2 < q < 00. The g-species annihilation process can 
be reformulated as a two-species annihilation model by 
combining g — 1 of the species into one group^_J_^) and 
the remaining specie into a second group (B) 30]. The 
reaction scheme becomes A + A Q and A + B — + 0. 
The ratio between the reaction rates of the two channels, 
, is a continuous parameter that need not necessarily 
correspond to an integer q. 

The transition rates arc as in the cyclic trapping re- 
action: (n, m,l, . . .) (n — 1, to — 1, /, . . .) occurs with 
rate nm/N. For symmetric initial conditions, the num- 
ber density p = a = b = c= -- - satisfies 



dp I -.s 2 



(23) 



and the initial condition p(0) = 1. The concentration is 
therefore 



[H-(q-l)t]-l. 



(24) 



Fluctuations can be obtained following the same 
straightforward steps the led to the evolution equations 
for the moments and we merely highlight the derivation. 
For q = 3, the probability distribution P(n, to, I) evolves 
according to ijHJl with the operator Cab defined via 



CabP = N-^ (A^Aij - 1) [nmP] 



(25) 



The probability distribution F{a, (3, 7) evolves according 
to (|ll|l with the evolution operator now being 

MabF = p{d^ + dp)[{a + (3)F] + ^{do.+ dpf F. (26) 

For arbitrary q, there are q{q — l)/2 such operators. 
Again, the first moment of F vanish; the second moments 
are coupled as follows 



dt 

d{aP) 
dt 



= -2{q - l)p{a^) - 2{q - l)p{af3) + {q - l)p\ 
= -2p{a^)-2{2q-?,)p{aP)+p\ (27) 



In contrast with the cyclic trapping reaction, the cross- 
correlation initially grows, although asymptotically it is 
again negative. 

Let U = (a2) + {q- l)(a/3) and V = (a^) - (a/3); 
the former quantity measures the correlation between n 
and the total particle number, the latter measures the 
correlation between n and n—m. To treat different values 
of q on the same footing, we rescale the time variable and 
introduce r = {q—l)t. The number density H23|) becomes 
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FIG. 2: The cross-correlation versus time for the cyclic trap- 
ping model and the three-species annihilation. 



p = (1 + r) ^ and the rate equations for the quantities 
U and V are 



dr 

ClT 9^1 



9-2 2 



(28) 



The solution for the first quantity is there- 
fore g-independent and apart from the numer- 
ical prefactor, as in the cyclic trapping model, 
U + r)-i - (1 + r)~4j rp^o different behaviors 

are found for the second quantity: 



(1 



(1 



£--2 

U(l + r)-iln(l+r) 



N-2 



q-2 



q = S 



(29) 



Asymptotically, the cross-correlation is negative be- 
cause (a/3) ~ "^^2^ ^^^^ ^'^ generically, fluctuations be- 
tween the numbers of different species are anti-correlated. 
Early on, the cross-correlation increases, but after a short 
transient it becomes negative (Fig. 2). 

In the long time limit, tr^ ~ (1 — q^^)V: 



3-q 



a 

P 



q<3; 
Int 9 = 3; 
0(1) q>3. 



(30) 



Therefore, fluctuations are relevant asymptotically only 
when g < 3. Applying the criterion y^Va^(t/) ^ Np{tf) 
yields the final time tf{N) and consequently, the typical 
final number of particles 



'iV(3-g)/2 g<3. 
0(1) q>3. 



(31) 



FIG. 3: The average total number of particles in the final 
state as a function of the system size for q — 4. The data 
represents an average over 10* realizations of the four-species 
annihilation process. 



There is an algebraic growth in the fluctuations domi- 
nated regime q < 3. At the critical point qc = 3, log- 
arithmic growth occurs. Otherwise, the final number of 
particles saturates at a finite value. Still, the final num- 
ber diverges, Nf ~ (g— 3)^^, in the vicinity of the critical 
point, q I 3. The saturation is illustrated in Fig. 3 using 
numerical simulation data for the g = 4 case. 

The case g = 2 is special since there is a conservation 
law (rt — m is conserved) and therefore, ^ = 0. Conse- 
quently, ^ p ^ t^^ . If a g-species aggregation, rather 
than annihilation process is considered, that is, when A 
and B interact the outcome is either A + B A ot 
A + B ^ B (both taken with equal probability), then 
this anomaly disappears [23, ^^"^ Ell- H31(l holds for 
g = 2 as well. 

One may ask "why is the critical number of species 
equal to 3?". Mathematically, the answer is ultimately 
related to the smaller eigenvalue of the 2x2 matrix cou- 
pling the second moments. Yet given the simplicity of 
the multispecies urn model, a heuristic and more illumi- 
nating derivation may be possible after all. Finding such 
an argument is an interesting challenge. 

In this context, we mention a generalization of the 
urn model from two-particle to the many-particle inter- 
actions. That is, instead of picking 2 particles, p particles 
are picked and if they are all of different species, all are 
removed from the system (this process is well defined for 
p < g) . The rate equations for the second moments yield 
the critical number of species (see Appendix) 



qdp) = 2p- 1. 



(32) 



Thus, for ternary interactions gc = 5. The structure of 
the phase transition remains the same. There is an alge- 
braic growth in the total number of remaining particles as 
a function of the system size when g < gc, a logarithmic 
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growth at the critical point q = qc, and saturation above 
the critical point q > qc- Last, we mention that a similar 
phase transition underlies two-species reaction processes 
of the type constructed from the g-species annihilation by 
separating species into two groups. In this case, although 
the transition depends on the reaction rates rather than 
the number of species, its structure remains the same. 

V. CONCLUSION 

In summary, we considered interacting particle systems 
undergoing depletion on the mean-field level. We showed 
that finite size fluctuations display a rich behavior. The 
behavior is rather generic and applies to a wide class of 
stochastic processes with a decreasing number of parti- 
cles. Typically, there is a phase transition as a function 
of the number of species or the reaction rates. In one 
region of parameter space, the final number of particles 
grows algebraically, and in the other, it saturates at a 
finite value. The critical case is marked by a logarithmic 
growth. We conclude that the final number of particles 
as a function of system size provides a practical probe of 
statistical fluctuations. 

The findings in the cyclic trapping model have a neat 
game theoretic implication. In a rock-paper-scissor game 
involving fixed strategy players and loser-is-out rules, the 
game ends when all remaining players have the same 



strategy. If players pair randomly, then the ultimate 
number of winners scales logarithmically with the total 
number of players. Intuitively, we expect that a similar 
law emerges for tournaments with simultaneous play. 

Several questions arise naturally. Can one explain 
the critical number of species using heuristic arguments? 
Statistical properties of the final state and how the sys- 
tem approaches it are interesting as well. For example, 
what is the number distribution of remaining particles? 
How different are statistical properties of the system at 
the very end of the process when only a single species re- 
mains? We observed that the convergence to the asymp- 
totic behavior is much faster when the first extinction 
occurs compared with the very end state when only a 
single species remains. Last, an interesting question in- 
volves extremal characterization 0, |23|: What is the 
probability that one of the species is always the most or 
least numerous? 

We studied systems undergoing depletion. However, 
there are processes in which depletion is possible but not 
certain, for example, infection processes It will be 

interesting to investigate finite size fluctuations in this 
related class of interacting particle systems. 

We are thankful to Matt Hasting, Sid Redner, and 
Zoltan Toroczkai for useful discussions. This research 
was supported by DOE(W-7405-ENG-36). 
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APPENDIX A: p-PARTICLE ANNIHILATION 



For p-particle annihilation, the density evolves according to 



(Al) 



The evolution operators for the Fokker-Planck equation for F are straightforward generalizations of Eq. H26|l. For 
example, for the ternary (p = 3) annihilation process A + B + C 0, 



MabcF = p\d^ + dp + d^)[{a + 13 + ^)F] + ti^{d^+dp + d.^ F. 



The second moments are coupled as follows 



dt 



o (^'-Dr^) VM>r^ la) 
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(A3) 



Introducing the time variable t — {p ~ the density is simply p ~ (1 + t) ^^'-p ^\ The quantity 

U = (a^) + {q — l){a(3) satisfies ^C/ + ^zjP^^^U — -^^p^ and the solution is again g-independent 



U{t) 



2p- 1 



(1 + r) p-i - (1 + t)" 



(A4) 



The quantity V = (a^) - {a(3) satisfies + (p-%^qli) P^~^V = (p-i)(g-i) P^- The solution reads 

F(r) 



(1 + t) p-i-(l + r) (p-i)(<!-i) 



q-(2p-l) 

5(^(1 + T)-i^ln(l + r) 



g = 2p - 1. 



(A5) 



Interestingly, in the fluctuation dominated regime, q < 2p — 1, the exponent governing the terminal time is p- 

g- 1 

independent, t/ ^ N^!~ . The final number of particles is 



N 2{p-i) q<2p-l; 
Nfr^^lnN q^2p-l; 

0{1) q>2p-l. 



(A6) 



